Down stream analysis of the dataset ‘4_A549_27plex’
knitr::opts_chunk$set(dev = "png",dpi = 80)
library(ggforce)
library(tidyverse)
dcolor <- ggsci::scale_color_d3('category20')
dfill <- ggsci::scale_fill_d3('category20')
ccolor <- viridis::scale_color_viridis()
style <- function(x,color,method='UMAP',axis=c(1,2),Theme=NULL,
axis.text=FALSE,axis.title=FALSE,coord_fix =TRUE,
palette=TRUE,title=NULL,
legend=TRUE,scale_color_log10=FALSE,
guide_pointSize=2.5,...) {
axis <- paste0(method,axis)
if(is.null(Theme)) Theme <- theme_bw
g <- ggplot() + Theme()
g <- g + geom_point(aes_(x=as.name(axis[1]),y=as.name(axis[2]),color=as.name(color)),x,...)
if(coord_fix) {
g <- g + coord_fixed()
}
if(!axis.text) {
g <- g + theme(axis.text = element_blank(),axis.ticks = element_blank())
}
if(!axis.title) {
g <- g + theme(axis.title = element_blank())
}
if(is.null(title)) {
g <- g + ggtitle(method)
}else if(!is.na(title)) {
g <- g + ggtitle(title)
}
if(!legend) g <- g + theme(legend.position = 'none')
if(scale_type(as.matrix(x[,color]))=='discrete' & !is.null(guide_pointSize)) {
g <- g + guides(color = guide_legend(override.aes = list(size=guide_pointSize)))
}
if(palette) {
if(scale_type(as.matrix(x[,color]))=='discrete') {
g <- g + ggsci::scale_color_d3('category20')
}else{
if(scale_color_log10){
g <- g + viridis::scale_color_viridis(trans='log10')
}else{
g <- g + viridis::scale_color_viridis()
}
}
}
return(g)
}
tib2df <- function(tib,RowNames=1) {
rn <- pull(tib,all_of(RowNames))
df <- tib %>%
select(-RowNames) %>%
as.data.frame()
rownames(df) <- rn
return(df)
}
data <- tibble(
path = list.files("4_A549_27plex/Raw/",full.names = T),
sampleName = c('0h-r1', '0h-r2', '0p5h-r1', '0p5h-r2', '24h-r1', '24h-r2', '48h-r1', '48h-r2')
) %>%
separate(sampleName,into=c('treat','rep'),remove = F) %>%
mutate(
metric = map(path,~{
.x %>%
read_csv(show_col_types=F) %>%
unite(cell,tileID,cellID,sep = '-',remove = F)
}),
IF = map(metric,~{.x[,-(2:8)]}),
meta = map2(metric,IF,~{
.x[,1:8] %>%
unite(cell,tileID,cellID,sep = '-',remove = F) %>%
mutate(totalExp = rowSums(.y[,-1]))
}),
) %>%
select(-path,-metric)
data %>%
select(-IF) %>% unnest(meta) %>%
ggplot(aes(x,y,color=treat)) +
theme_void() +
geom_point(size=.1) + dcolor +
facet_grid(rep~treat) + coord_fixed()
data %>%
select(sampleName,meta) %>% unnest(meta) %>%
ggplot(aes(sampleName,totalExp)) + geom_violin() + geom_boxplot()
dataf <- data %>%
mutate(
idx = map(meta,~{.x$totalExp > 1.3}),
IF = map2(IF,idx,~{.x[.y,]}),
meta = map2(meta,idx,~{.x[.y,]})
) %>%
mutate(
idx2 = map(IF,~{rowSums(.x[,-1]==0) == 0}),
IF = map2(IF,idx2,~{.x[.y,]}),
meta = map2(meta,idx2,~{.x[.y,]}),
nCell = map(idx2,sum) %>% unlist()
) %>%
select(-idx,-idx2)
dataf %>%
select(sampleName,meta) %>% unnest(meta) %>%
ggplot(aes(sampleName,totalExp)) +
geom_violin() + geom_boxplot()
dataf %>%
ggplot(aes(sampleName,nCell)) +
theme_classic() +
geom_bar(stat='identity') +
scale_x_discrete(limits=rev) +
theme(axis.title.y = element_blank()) +
coord_flip() + theme(aspect.ratio = .35)
tmp <- dataf %>% select(sampleName,IF) %>% unnest(IF)
idx <- with(tmp,model.matrix(~0+sampleName))
colnames(idx) <- sub("sampleName","",colnames(idx))
avgmat <- tmp[,-(1:2)] %>% as.matrix() %>% t %>% {t(.%*%idx)/colSums(idx)} %>% scale %>% asinh %>% t
pheatmap::pheatmap(avgmat,cluster_cols = F,cellwidth = 10,color=viridis::viridis(10),show_rownames = T)
Meta <- dataf %>%
select(sampleName,treat,rep,meta) %>% unnest(meta) %>%
unite(Cell,treat,rep,cell,sep = '-',remove = F)
X <- bind_rows(dataf$IF) %>% select(-cell) %>%
scale() %>% asinh()
rownames(X) <- Meta$Cell
p <- prcomp(X)
screeplot(p,type='l',n=30)
total <- bind_rows(dataf$meta) %>% pull(totalExp)
cor(p$x,total) %>%
plot(xlab = 'PC',ylab = 'CorCoef_vsTotalExp')
um <- uwot::umap(p$x[,2:6],
metric = 'cosine',
n_neighbors = 30,
min_dist = .1) %>%
as_tibble(rownames = 'Cell') %>%
dplyr::rename(UMAP1=2,UMAP2=3)
Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if `.name_repair` is omitted as of tibble 2.0.0.
Using compatibility `.name_repair`.
Metaum <- inner_join(Meta,um,by = 'Cell')
style(Metaum,'treat',title = NULL,size = .1,alpha = .6)
Warning: `aes_()` was deprecated in ggplot2 3.0.0.
Please use tidy evaluation idioms with `aes()`
style(Metaum,'sampleName',title = NA,size = .1,alpha = .6)
style(Metaum,'totalExp',title = NA,size = .1,alpha = .6)
Using the datasets of 0h and 0.5h treatment
Meta2 <- Meta %>% filter(treat %in% c('0h','0p5h'))
X2 <- dataf %>%
filter(treat %in% c('0h','0p5h')) %>%
pull(IF) %>% bind_rows() %>%
select(-cell) %>%
scale() %>% asinh()
rownames(X2) <- Meta2$Cell
p2 <- prcomp(X2)
screeplot(p2,type='l',n=30)
total2 <- Meta2$totalExp
cor(p2$x,total2) %>%
plot(xlab = 'PC',ylab = 'CorCoef_vsTotalExp')
um2 <- uwot::umap(p2$x[,2:6],
metric = 'cosine',
n_neighbors = 30,
min_dist = .1) %>%
as_tibble(rownames = 'Cell') %>%
dplyr::rename(UMAP1=2,UMAP2=3)
Metaum2 <- inner_join(Meta2,um2,by = 'Cell')
style(Metaum2,'treat',title = NULL,size = .1,alpha = .6)
ph2 <- phateR::phate(p2$x[,2:6],
mds.dist.method = 'cosine',
knn = 20)
Calculating PHATE...
Running PHATE on 19597 observations and 5 variables.
Calculating graph and diffusion operator...
Calculating KNN search...
Calculated KNN search in 2.02 seconds.
Calculating affinities...
Calculated affinities in 0.10 seconds.
Calculated graph and diffusion operator in 2.14 seconds.
Calculating landmark operator...
Calculating SVD...
Calculated SVD in 1.36 seconds.
Calculating KMeans...
Calculated KMeans in 2.33 seconds.
Calculated landmark operator in 4.20 seconds.
Calculating optimal t...
Automatically selected t = 26
Calculated optimal t in 2.87 seconds.
Calculating diffusion potential...
Calculated diffusion potential in 0.58 seconds.
Calculating metric MDS...
Calculated metric MDS in 5.23 seconds.
Calculated PHATE in 15.02 seconds.
Metaum2 <- Metaum2 %>%
mutate(PHATE1 = ph2$embedding[,1],PHATE2 = ph2$embedding[,2])
style(Metaum2,'treat','PHATE',size=.1)
style(Metaum2,'rep','PHATE',size=.1,palette = F)
set.seed(1); km <- kmeans(ph2$embedding,centers=5)
Metaum2 %>%
mutate(km=as.factor(km$cluster)) %>%
style('km','PHATE',size=.1)
library(slingshot)
Loading required package: princurve
Loading required package: TrajectoryUtils
Loading required package: SingleCellExperiment
Loading required package: SummarizedExperiment
Loading required package: MatrixGenerics
Loading required package: matrixStats
Attaching package: ‘matrixStats’
The following object is masked _by_ ‘.GlobalEnv’:
count
The following object is masked from ‘package:dplyr’:
count
Attaching package: ‘MatrixGenerics’
The following objects are masked from ‘package:matrixStats’:
colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse, colCounts, colCummaxs, colCummins, colCumprods, colCumsums, colDiffs, colIQRDiffs, colIQRs, colLogSumExps,
colMadDiffs, colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats, colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds, colSums2, colTabulates,
colVarDiffs, colVars, colWeightedMads, colWeightedMeans, colWeightedMedians, colWeightedSds, colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet, rowCollapse,
rowCounts, rowCummaxs, rowCummins, rowCumprods, rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps, rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks, rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars, rowWeightedMads, rowWeightedMeans,
rowWeightedMedians, rowWeightedSds, rowWeightedVars
Loading required package: Biobase
Welcome to Bioconductor
Vignettes contain introductory material; view with 'browseVignettes()'. To cite Bioconductor, see 'citation("Biobase")', and for packages 'citation("pkgname")'.
Attaching package: ‘Biobase’
The following object is masked from ‘package:MatrixGenerics’:
rowMedians
The following objects are masked from ‘package:matrixStats’:
anyMissing, rowMedians
Attaching package: ‘SummarizedExperiment’
The following object is masked from ‘package:Seurat’:
Assays
The following object is masked from ‘package:SeuratObject’:
Assays
sce <- getLineages(ph2$embedding,as.character(km$cluster),start.clus = '5') %>%
getCurves() %>% slingPseudotime() %>%
as_tibble(rownames='Cell') %>% mutate(Lineage1 = Lineage1/max(Lineage1)) %>%
rename(PseudoTime = Lineage1)
Warning: useNames = NA is deprecated. Instead, specify either useNames = TRUE or useNames = TRUE.Warning: useNames = NA is deprecated. Instead, specify either useNames = TRUE or useNames = TRUE.Warning: useNames = NA is deprecated. Instead, specify either useNames = TRUE or useNames = TRUE.Warning: useNames = NA is deprecated. Instead, specify either useNames = TRUE or useNames = TRUE.Warning: useNames = NA is deprecated. Instead, specify either useNames = TRUE or useNames = TRUE.
Metaum2 <- left_join(Metaum2,sce,by='Cell')
Metaum2 %>% style("PseudoTime","PHATE",size=.2) + theme_void()
Metaum2 %>% style("PseudoTime","UMAP",size=.2) + theme_void()
Metaum2 %>%
ggplot(aes(PseudoTime,fill=treat)) + geom_histogram(position = 'dodge') + theme_classic() + dfill + coord_fixed(ratio=.0003) +
theme(axis.ticks = element_blank(),panel.grid = element_blank())
LOESS <- dataf %>%
filter(treat %in% c("0h","0p5h")) %>%
select(sampleName,IF) %>% unnest(IF) %>%
unite(Cell,sampleName,cell,sep='-') %>%
left_join(select(Metaum2,Cell,PseudoTime),.) %>%
gather(Prot,ExpLevel,-(1:2)) %>%
group_by(Prot) %>% nest() %>%
mutate(LOESS = map(data,~loess(ExpLevel~PseudoTime,.x,span = 0.2)),
PRED = map(LOESS,~predict(.x,newdata=seq(0,1,by=0.01))),
newdata = map(PRED,~tibble(PseudoTime=seq(0,1,by=0.01),PRED=.x)) )
Joining with `by = join_by(Cell)`
loessPl <- function(x,ratio = .35){
ggplot(x,aes(PseudoTime,scaled,color=Prot)) +
geom_line() + theme_bw() +
theme(legend.position = 'none') +
theme(axis.ticks.y = element_blank(),
strip.background = element_rect(fill=NA,color=NA),
panel.grid.major.x = element_blank(),panel.grid.minor = element_blank()) +
coord_fixed(ratio=ratio) + labs(y = 'scaled expression level')
}
LOESS_scaled <- LOESS %>% select(Prot,newdata) %>%
mutate(newdata = map(newdata,~mutate(.x,scaled = (PRED-min(PRED))/(max(PRED)-min(PRED)) ))) %>%
unnest(newdata) %>% ungroup()
usep <- c("p4EBP1/2/3","pATM","pELF2A","pERK","pFAK","pp38","pSMAD1/5","pSMAD3","pSTAT3")
LOESS_scaled %>% filter(Prot %in% usep) %>% loessPl()
LOESS_scaled %>% filter(Prot %in% usep) %>% loessPl() + facet_wrap(~Prot)
LOESS_scaled %>% filter(!Prot %in% usep) %>% loessPl() + facet_wrap(~Prot)
LOESS_scaled %>% filter(!Prot %in% usep) %>% loessPl()
sessionInfo()
R version 4.2.2 (2022-10-31)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.2.1
Matrix products: default
LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods base
other attached packages:
[1] slingshot_2.6.0 TrajectoryUtils_1.6.0 SingleCellExperiment_1.20.1 SummarizedExperiment_1.28.0 Biobase_2.58.0 MatrixGenerics_1.10.0 matrixStats_1.0.0
[8] princurve_2.1.6 lubridate_1.9.2 forcats_1.0.0 stringr_1.5.0 dplyr_1.1.2 purrr_1.0.1 readr_2.1.4
[15] tidyr_1.3.0 tibble_3.2.1 tidyverse_2.0.0 ggforce_0.4.1 ggplot2_3.4.2 Seurat_4.3.0.1 SeuratObject_4.1.3
[22] sp_2.0-0 GenomicRanges_1.50.2 GenomeInfoDb_1.34.9 IRanges_2.32.0 S4Vectors_0.36.2 BiocGenerics_0.44.0
loaded via a namespace (and not attached):
[1] utf8_1.2.3 spatstat.explore_3.2-1 reticulate_1.30 tidyselect_1.2.0 RSQLite_2.3.1 AnnotationDbi_1.60.2 htmlwidgets_1.6.2
[8] grid_4.2.2 Rtsne_0.16 munsell_0.5.0 codetools_0.2-19 ica_1.0-3 future_1.33.0 miniUI_0.1.1.1
[15] withr_2.5.0 spatstat.random_3.1-5 colorspace_2.1-0 progressr_0.13.0 knitr_1.43 rstudioapi_0.15.0 ROCR_1.0-11
[22] tensor_1.5 listenv_0.9.0 labeling_0.4.2 GenomeInfoDbData_1.2.9 polyclip_1.10-4 bit64_4.0.5 farver_2.1.1
[29] pheatmap_1.0.12 rprojroot_2.0.3 parallelly_1.36.0 vctrs_0.6.3 generics_0.1.3 xfun_0.39 timechange_0.2.0
[36] R6_2.5.1 phateR_1.0.7 bitops_1.0-7 spatstat.utils_3.0-3 cachem_1.0.8 DelayedArray_0.24.0 promises_1.2.0.1
[43] scales_1.2.1 vroom_1.6.3 gtable_0.3.3 globals_0.16.2 goftest_1.2-3 rlang_1.1.1 pamr_1.56.1
[50] splines_4.2.2 lazyeval_0.2.2 spatstat.geom_3.2-4 yaml_2.3.7 reshape2_1.4.4 abind_1.4-5 stylo_0.7.4
[57] httpuv_1.6.11 tools_4.2.2 tcltk_4.2.2 ellipsis_0.3.2 jquerylib_0.1.4 RColorBrewer_1.1-3 ggdendro_0.1.23
[64] proxy_0.4-27 ggridges_0.5.4 Rcpp_1.0.11 plyr_1.8.8 sparseMatrixStats_1.10.0 zlibbioc_1.44.0 RCurl_1.98-1.12
[71] deldir_1.0-9 pbapply_1.7-2 viridis_0.6.4 cowplot_1.1.1 zoo_1.8-12 ggrepel_0.9.3 cluster_2.1.4
[78] here_1.0.1 magrittr_2.0.3 data.table_1.14.8 scattermore_1.2 lmtest_0.9-40 RANN_2.6.1 fitdistrplus_1.1-11
[85] hms_1.1.3 patchwork_1.1.2 mime_0.12 evaluate_0.21 xtable_1.8-4 gridExtra_2.3 compiler_4.2.2
[92] KernSmooth_2.23-22 crayon_1.5.2 htmltools_0.5.5 later_1.3.1 tzdb_0.4.0 DBI_1.1.3 tweenr_2.0.2
[99] MASS_7.3-60 Matrix_1.5-4.1 cli_3.6.1 parallel_4.2.2 igraph_1.5.0.1 pkgconfig_2.0.3 plotly_4.10.2
[106] spatstat.sparse_3.0-2 bslib_0.5.0 XVector_0.38.0 digest_0.6.33 sctransform_0.3.5 RcppAnnoy_0.0.21 tsne_0.1-3.1
[113] spatstat.data_3.0-1 Biostrings_2.66.0 rmarkdown_2.23 leiden_0.4.3 uwot_0.1.16 DelayedMatrixStats_1.20.0 shiny_1.7.4.1
[120] lifecycle_1.0.3 nlme_3.1-162 jsonlite_1.8.7 viridisLite_0.4.2 fansi_1.0.4 pillar_1.9.0 ggsci_3.0.0
[127] lattice_0.21-8 KEGGREST_1.38.0 fastmap_1.1.1 httr_1.4.6 survival_3.5-5 glue_1.6.2 png_0.1-8
[134] bit_4.0.5 tcltk2_1.2-11 class_7.3-22 stringi_1.7.12 sass_0.4.7 blob_1.2.4 memoise_2.0.1
[141] irlba_2.3.5.1 e1071_1.7-13 future.apply_1.11.0 ape_5.7-1